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I.  INTRODUCTION 


We  consider  the  problem  of  determining  the  probability  density  func¬ 
tion  of  the  output  of  an  RC  filter  when  the  input  is  a  binary  random  process. 

Such  problems  have  been  around  for  forty  or  so  years,  and  arose  during  that 
time  largely  out  of  intrinsic  interest  in  the  development  of  the  theory  of 
random  processes.  Their  solutions  and  the  methods  used  to  solve  them  have 
provided  insight,  and  have  led  to  many  useful  applications  [l,  2]  .  More 
recently,  filtered  binary  processes  have  received  attention  as  theoretical 
models  of  a  single  channel  in  a  nerve  membrane  distorted  by  the  low- pass 
filtering  action  of  recording  equipment  [5,  15]  .  Most  previous  work  has  focused 
attention  on  the  special  case  in  which  the  binary  input  process  is  the  random 
telegraph  signal  (with  exponentially  distributed  intervals),  and  for  this  case 
results  are  available  for  the  output  density  function  as  well  as  for  related 
quantities  such  as  the  interval  statistics  and  level  crossings  of  the  output 
[5,9-11,  13,  14,  17].  Few  results  have  been  obtained  in  the  more  difficult  case 
in  which  the  intervals  of  the  input  are  not  exponentially  distributed. 

For  independent  and  identically  distributed  intervals  with  arbitrary 
statistics,  McFadden  [6]  derived  integral  equations  from  which  the  output 
density  can  be  obtained.  However,  because  of  the  complexity  of  these  integral 
equations,  he  was  able  to  find  a  solution  in  only  one  special,  non-trivial  case. 
Another  case  that  has  received  considerable  attention,  but  with  limited  success, 
is  that  in  which  the  binary  process  is  the  result  of  hard-limiting  of  a  stationary 
Gaussian  noise  with  exponential  correlation  [4,7,9].  In  this  case,  the 
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intervals  of  the  binary  process  are  no  longer  statistically  independent.  Al¬ 
though  exact  results  are  available  for  only  one  special  ratio  of  time  constants 
[4],  approximations  have  been  obtained  by  use  of  the  Fokker-Planck  equa¬ 
tion  [9]  . 

In  this  paper,  attention  will  be  confined  to  binary  processes  of  the 
type  considered  by  McFadden;  i.e.,  each  binary  process  will  be  assumed  to 
be  constructed  from  an  equilibrium  renewal  process  [3,  16]  so  that  its  intervals 
will  be  independent  of  one  another  and  of  the  state  of  the  binary  process  at  the 
transition  times.  We  will  extend  McFadden's  work  in  several  directions: 

(i)  by  generalizing  the  binary  input  process  to  have  different  probability 
density  functions  for  the  up  and  down  intervals,  (ii)  by  finding  approximations 
based  upon  the  Fokker-Planck  equation,  fiii )  by  deriving  new  integral  equa¬ 
tions  for  the  relevant  density  functions,  (iv)  by  developing  matrix  and  iterative 
methods  for  solutions  of  the  integral  equations,  and  (v)  by  investigating  trans¬ 
formations  of  the  integral  equations  into  differential  equations.  Some  numeri¬ 
cal  results  which  compare  the  various  approaches  will  be  given. 

The  paper  begins  in  the  next  section  with  a  brief  review  of  the  system 
model,  the  integral  equations  of  McFadden,  and  a  summary  of  known  results. 
Sec.  Ill  discusses  the  moments  of  the  filtered  process,  and  Sec.  IV  gives 
approximations  based  upon  the  Fokker-Planck  equation.  In  Sec.  V,  new 
integral  equations  are  derived,  and  methods  for  their  solution  considered. 
Transformation  of  the  integral  equations  into  differential  equations  is 
examined  in  Sec.  VI,  numerical  results  are  presented  in  Sec,  VII,  and  the 
final  section  summarizes  and  discusses  the  results. 
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II.  SYSTEM  MODEL  &  KNOWN  RESULTS 


This  section  gives  a  brief  summary  of  the  system  model  and  known 
results,  and  also  serves  to  define  the  various  quantities  that  enter  into  the 
analyses. 

The  system  is  governed  by  the  differential  equation 

-  +  3 y(t)  =  3x(t)  ,  3>0  fl) 

where  3  ^  is  the  RC  filter  time  constant,  x(t)  =  i  1  is  the  binary  input 
process  and  y(t),  |  y (t )  |  <  1,  is  the  filter  output.  We  shall  frequently  use 


A.  The  Binary  Input  Process.  The  binary  process  x(t)  is  characterized  by 
the  time  intervals  between  its  transitions,  which  are  taken  to  be  independent 
random  variables.  Further,  the  time  intervals  corresponding  to  the  x(t)  =  -1 
state  will  be  taken  to  be  identically  distributed  with  one  probability  density 
function,  while  those  corresponding  to  the  x(t)  =  *1  state  to  be  identically 
distributed  with  a  different  density.  These  densities  and  some  quantities  later 
needed  are 


f q (t )  =  p.d.  f.  of  x(t)  =  -1 

interval  length 

(2a) 

f  j  (t)  =  p.  d.  f.  of  x(t)  =  *1 

interval  length 

(2b) 

X 

/- 

^  .  -  /  t  f  .  (t )  dt 

(2c ) 

1  o  1 

X 

F . ( s )  =  £  Cf  . (t ) }  =  /  e’St 

*1  IV 

f  .  (t )  dt 

(2d) 
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1 
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5  .  (t )  =  /  f  .  (t ' )  dt ' 
i 

(2e) 
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zero 


All  of  the  f's  are  to  be  considered  as  being  defined  on  [0,  *  ] ,  and  to  be 
for  negative  arguments,  U.  is  the  mean  value  of  the  up-  or  down  interval, 
FJs)  the  Laplace  transform  of  the  interval  density,  and  ~.(t)  the  comple¬ 
ment  of  the  interval  probability  distribution  function.  The  quantities  in 
(2a)- (2e)  will  sometimes  be  written  without  subscripts  in  the  symmetric  case 
in  which  f(J(t)-f1(t)-fCt),/*-0-/*1-/i,  F0(s)=F1(s)=F(s) ,  and  U)=J(t)  . 

B.  The  Output  Process.  The  output  process  y(t)  consists  of  segments  of 
rising  and  decaying  exponentials,  and  consequently  will  have  local  minima  and 
maxima  at  the  transition  points  of  x(t).  The  probability  density  functions  of 
y(t)  at  these  transition  points  enter  into  the  analyses.  At  this  point,  it  is 
convenient  to  define  six  different  probability  density  functions  associated  with 
y(t),  and  these  are 

pQ(y)=  p.d.f.  of  y(t)  at  a  minimum  point  (3a) 

p ^  (y )  =  p.d.f.  of  y(t)  at  a  maximum  point  (3b) 

p  (y)  =  p.d.f.  of  y(t )  at  a  time  picked  at  random  during  an 

x(t)  =  -1  interval  (3c) 

pjy)  =  p.d.f.  of  y(t)  at  a  time  picked  at  random  during  an 

x(t)  =  +1  interval  (3d) 

p(y)  =  unencumbered  p.d.f.  of  y(t)  at  a  time  picked  at  random  (3e) 


p j,p(y)  =  Fokker-Planck  approximation  to  p(y) 


All  of  the  p's  are  zero  outside  the  interval  [-1,1]  -  this  will  be  tacitly 
assumed  throughout  and  will  not  be  stated  each  time  an  expression  for  one  of 
the  p'3  is  given.  In  the  symmetric  case,  we  have  Pg(y)  =  p^(-y)  and 

p  (y)  =  ?J-y). 


Some  elementary  relations  between  p  (y),  p^(y)  and  p(y)  are 


p(y) 

P  (y ) 

p_(y> 
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-A. 

“0 

“l 

4* 

0 

“l 

2- 

0 

U  4“ 

U 

o 

l 

?_(y)  -  — 4— 
0  1 


(1  -  y)  p(y) 


(1  -  y)  p(y) 


?Jy) 


(4a) 


(4b) 


(4c) 


Equations  (4b)  and  (4c)  follow  from  a  theorem  on  conditional  expectation  of 
Mazo  and  Salz  [8]  which  says  that  E  [  y  j  y]  =  0.  Applying  this  to  (1)  gives 
E  [x  (t )  !  y  (t )  ]  =  y(t )  which, in  turn.implies  (4b)  and  (4c)  (see  the  derivation  in 
[ll]  for  a  similar  situation).  From  (4b)  and  (4c)  it  follows  that 

u0(l+y)p_(y)  =  u1(l-y)p+(y)  (4d) 

which  also  can  be  obtained  by  equating  the  average  number  per  second  of  the 
upward  and  downward  crossings  of  the  level  y. 

C.  The  Basic  Integral  Equations.  The  basic  integral  equations  relating  the 
interval  densities  and  densities  of  the  output  process  are.  with  T  =  1/3  , 


P0(y>=  rh  /  d-p1Mb(T!nr 


+  n 


(5a) 


?i(y)  =  —  j  d"  Po(")£ilTta  ~ 


1  -_2 

y 


p>> 1  /  d"  pi(")  5o(T  M  Hr 


u  a  -y) 

0  y 


p+<yl  ~  /  j  d'p0(‘)3i(T  5,1  Hr 


(5b ) 


(5c ) 


(5d) 


The  precise  mathematical  conditions  required  are  satisfied  for  the  processes 
considered  here. 


Equations  (5b)  and  l5d)  were  derived  by  McFadden  (cf.  (5)  and  (7)  of  [6]),  who 
wrote  them  in  a  slightly  different  form,  and  are  the  result  of  considering  the 
output  y(t)  over  a  time  interval  during  which  x(t)  =  *1.  Exactly  following  the 
same  steps  for  an  x(t)  =  -1  interval  leads  to  (5a)  and  (5c).  In  the  symmetric 
case,  McFadden  changed  the  independent  variable  from  y  to  t  by  means  of 


v  =  1  -  Ze 


and,  by  Laplace  transforming  with  respect  to  the  new  variable  t,  transformed 
(5b)  and  (5d)  into  algebraic  expressions  which  enabled  him  to  obtain  a  solution 
in  one  special  case.  To  do  the  same  type  of  transformation  in  connection  with 
(5a)  and  (5c)  would  require  the  slightly  different  change  of  variable 


y  =  -1  ^  2e 


D.  Known  Results.  The  only  nontrivial  cases  in  which  exact  results  are  known 
are  those  of  exponentially  distributed  intervals  and  McFadden's  special  case. 

A  brief  summary  of  these  follows. 

(i)  Symmetric  Case:  f(t)  =  ae  at  ,  [10,  17] 


Z2a  B (a.  a) 


=  aT 


p  (y)  =  (1  -  y)  p (y ) 


where  B(u,  v)  =  7(u)  7(v)/T(u  +  v)  is  the  beta  function. 


( 9a ) 


HI,  MOMENTS 


The  momenta  E  [yn(t)]  ,  n  =  1 ,  2,  .  .  .  can  often  be  evaluated  in  situa¬ 
tions  of  the  present  kind  without  first  getting  the  density  p(y)  [7 ]  .  In  the 
symmetric  case,  all  of  the  odd  moments  are  zero.  In  general,  the  moments 
can  be  expressed  in  terms  of  the  Laplace  transforms  of  the  interval  densities, 
F^(s)  and  F^(s),  and  the  first  two  nonzero  moments  in  the  symmetric  and 
asymmetric  cases  are: 


(i)  Symmetric  Case: 


2  .  2  1  -  F(S) 

y  =  1  - 


3u  1  +  F(S) 


33u  [l  +F(=)][l  -  F(25)][l  rF(33)] 


(ii)  Asymmetric  Case: 


U1  -uo 


u  T  u 
0  1 


1  - 


[1  -  Fq(3  )][  1  -  F  L  (3  )] 


3(ao  +  ur 


1  -Fo(S)Fi(?) 


(11a) 


4  ,  8  [l  -  2F(3)  +  2F(23)  -  F(-)F(23)][l  -  F(33)] 

y  =  1  -  —  - 7 - —  :-~-r - : _ t  r - c— — — - (1  lo) 


(11c) 


(lid 


Equation  (11a)  is  derived  in  the  Appendix,  and  its  generalization  to  (lid)  is 
also  outlined  there,  A  recursive  method  for  obtaining  all  of  the  moments  of 
y(t)  as  well  as  all  of  the  conditional  moments  E[yn(t)ix(t)  =  =  l1,  has  recentl 
been  developed  by  A.  Munford  [l9],  and  (11b)  has  been  deduced  from  his 
work.  The  n-th  moment  is  a  function  of  F^(k3)  and  F^(k3)  for 
k  =  1 ,  2,  .  .  .  n-  1 ;  and  becomes  increasingly  complicated  as  the  order  of  the 
moment  increases.  The  moments  are  significantly  more  complicated  in  the 
asymmetric  case  than  in  the  symmetric  case. 


IV. 


FCKKER-  PLANCK  A P PRC X IALA. T IC NS 


Apor  oximations  to  p(y)  based  upon  a  certain  linearity  assumption  can 
be  obtained  by  use  of  the  Fokker  -  Planck  equation.  The  Fokker- Planck  approx 
imations  are,  in  fact,  exact  in  the  case  that  x(t)  has  exponentially  distributed 
intervals  [5,  10],  and  are  close  approximations  when  x(t)  is  hard-limited 
Gaussian  RC  noise  [T ,  9]  .  It  will  later  be  shown  that  the  Fokker  -  Planck 
approximation  is  also  exact  in  the  case  of  the  McFadden  interval  density  with 
b  =  a  +  1  and  b  =  a  +  2.  In  some  other  cases,  the  approximations  are  accu¬ 
rate  to  within  a  few  percent.  Also,  the  approximations  provide  starting  points 
for  iterative  solutions  to  the  integral  equations,  to  be  considered  later.  Our 
use  of  the  Fokker  -  Planck  equation  closely  parallels  that  in  [9]  . 


A.  The  Symmetric  Case.  The  (extended)  Fokker  -  Planck  equation  for  p(y)  is 
readily  shown  to  be  [Q,  11] 


_d_ 

dv 


1,2  2 

—  ~  fl  -  y  )  ply) 


-  A(y)  p  (y)  =  0 


(12) 


in  which  A(y)  is  defined  as  the  limit  of  a  conditional  expectation  as 


A(y )  =  lim  — —  E  [y(t  ~  A )  -  y  ft )  y<t )] 
A  -»  0  A“ 


(13) 


In  general,  it  is  not  known  how  to  evaluate  A(y).  An  assumption  that  has  led 
to  an  excellent  approximation  in  one  case  [Q]  is  that  A(y)  is  approximately  a 
linear  function  of  v,  and,  for  that  reason,  we  make  the  same  assumption  here 


A(y )  = 


constant 


’  1  ■*) 


c  uy  ,  u  - 


in  (111  and  solving  the  ensuing  equation  then  gives 


~>  u- 


pFP(y'  = 


(1  -  V"> 


( 1: 


2u-  1  „ 

2  B(u,u) 


The  subscripts  'FP"  on  the  density  function  are  intended  to  emphasize  the 
fact  that  the  density  function  is  an  approximation  based  upon  the  Fokker- 

Planck  equation,  and  not  necessarily  an  exact  result.  The  unknown  constant 

~ 

u  can  be  found  by  evaluating  y“  by  two  different  ways  (see  (11a)  and  the 
Appendix),  and  comes  out  to  be 

1  -  F  (3) 


u 


Su[l  +F(6)]  -  2  [1  -  F  (S )  ] 


'  1  o ) 


It  will  later  be  shown  in  Sec.  VI-A  (see  (43)  et  seq.  )  that  an  exact 


expression  for  A(y)  in  terms  of  p(y)  and  p  (y)  is 


Alyl  .  Vy‘  ~  P0f-y’ 


2c u  d (y) 


(IT) 


which  provides  a  way  of  comparing  the  assumption  (14)  with  the  exact  A(y) 
when  p(y)  and  p^fy)  are  known. 

3.  The  Asymmetric  Case.  The  Fokker  -  Planck  approximation  in  the  asym¬ 
metric  case  follows  in  much  the  same  way  as  in  the  symmetric  case,  and  so 


the  details  will  not  be  given.  The  results  are 


w  he  r  e 
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LI, 


U 


(19) 


u  = 


1 

r 


u  ~  v  x  (u  -  v) 
u  v  ~  (u  ~  V ) 


(20) 


(21) 


and,  the  analog  of  (17)  is 

A(y)  =  PQ(y)  ~  <Vyl 

a  2  3  (UQ  A  u  L  )  p(y) 


(22) 


Equation  (lid)  is  an  alternate  expression  for  y  which  can  be  used  in  (20)  to 
give  the  value  of  the  parameter  u  in  terms  of  the  Laplace  transforms  of  the 


interval  densities. 


V.  INTEGRAL  EQUATION  APPROACH 


A.  Some  New  Integral  Equations.  In  this  section,  we  explore  ways  of  further 
using  the  basic  integral  equations  (5a)  -  (5d).  Equations  (5a)  and  (5b)  can  be 
used  in  the  integrands  of  one  another  to  give  integral  equations  in  either  p  (y) 
or  p.(y)  only;  viz.  . 


By  interchanging  orders  of  integration,  these  can  also  be  written  as 


where  the  kernels  are 


K  i  v 
0Y> 


1  - 


p0(y)  =  J  K0(y,")  p0(^)d’ 


Pl  (y)  =  J  Kj  (y,  ")  pj  (")  dv 


1 


fo(T  :nT^ 


Kiy,~)  = 


1  -  y 


maxiy,  ’") 
min(y,  “) 


■olT  ITT 


fi*T  * 


- 1 


(24 


(24 


(25 


(25 


Eqs.  (24a)  and  (24b  *  were  derived  independently  by  R.  FitzHugh  in  a  private 
correspondence. 


The  forms  of  (24a)  and  (24b)  lend  themselves  more  directly  to  numerical  tech¬ 
niques  than  do  (23a)  and  (23b).  A  matrix  technique  for  solving  them  will  be 
discussed  later  in  this  section. 


A  second  way  of  using  the  set  of  basic  integral  equations  is  that  of 
inverting  (5c)  and  (5d)  to  express  p^(y)  and  p^(y)  as  integrals  of  p  (y)  and 

p+(y): 


(26a) 


in  which 


Pj 


T 


y 


d~  p  ( ^ ) 


'-1 


(26b) 


h.  (t)  =  £ 


|  »ri(.)  j 

I  1-Fi(»)l 


(27) 


Equation  (2ba)  can  be  obtained  by  using  the  change  of  variable  (7)  in  (5a)  and 
(5c),  and  then  Laplace  transforming  both  with  respect  to  the  new  variable.  The 
resulting  two  equations  can  be  solved  for  the  transform  of  Pg(y)  in  terms  of  the 
transform  of  p  (y).  Inver s ion  then  leads  to  (26a).  Equation  (26b)  can  be  ob¬ 
tained  in  the  same  way  by  starting  with  (5b)  and  (5d)  and  using  the  change  of 
variable  (6). 

As  with  f  (t),  whenever  h^(t)  is  used,  it  will  be  tacitly  assumed  to  be 
zero  for  negative  t.  When  (26a)  and  (26b)  are  substituted  into  the  integrands  of 
(5c)  and  (5d),  the  results  are  two  equations  containing  only  p  (y)  and  p  (y). 
Employing  (4b)  and  (4c),  the  first  of  these  gives  an  equation  involving  only  the 
unencumbered  densitv  o(v);  i.e., 


p(y)  = 


— i 7  T  v/i 

l-y2  ,l+?  '-y  A 

-  1  -3 


d-  n-r)T,(-)  h  It  zn 


i  +  3 


(2 


which  can  also  be  written  as 


X 

p(y)  =  j  K(y,  "(pM  d’ 


(2 


where 


^  min(y,’1) 

K(y, - 1  -  — j—  f  "0(T&^-)  <3 

■'-l 


1  - 


Any  of  (23a),  23b)  or  (29)  can  be  solved  by  the  matrix  method  of  the  next  sub¬ 
section. 

The  h(t)  functions  will  usually  be  more  complicated  than  the  interval 
densities,  however,  there  are  some  cases  in  which  the  h(t)  are  particularly 
simple.  The  Laplace  transform  of  the  McFadden  interval  density  (10a)  is 


F(s) 


r  fo)  r  (sT  +  a) 

r  (a)  r  (sT  +  b) 


(3 


and  the  Laplace  transform  of  the  gamma  interval  density 


f(t) 


n+1  n  -at 
a  t  e 


n  =  0.  1. 


is 


F(s )  = 


(3 


Cases  in  which  these  lead  to  simple  forms  for  h(t)  are: 


(i)  McFadden,  b  =  a  +  1 


h(t)  =  a3  6(t) 


(34a) 


(ii)  McFadden,  b  =  a +  2 


h(t)  =  a3(a  +  l)  e'(2arl)3t  (34b) 

(iii)  McFadden,  b  =  a  +  3 

,  2a3  (a  +  l)(a+2)  -3(a+l)3t/2  .  ,  j~  2  ^  tr,  i  \ 

h(t)  =  —  !  —  — . -  e  sin  (V  3a  +oa  -  1  3t/Z)  (34c ) 

J . 3a ^  +  6a  -  1 

(iv)  Gamma,  n  =  0 

h(t)  =  a  5  (t)  (34d) 


(v)  Gamma,  n  =  1 

h(t)  =  a^  e  (34e) 

(vi)  Gamma,  n  =  2 

h(t)  =  {2  a?  /  */~3)  e  sin(VTat/2)  (34f) 

For  either  the  McFadden  density  with  b  -  a  an  integer,  or  the  gamma  density, 
l  h(t):  is,  in  general,  the  ratio  of  two  polynomials  in  s.  An  explicit  form 
for  h(t)  as  the  sum  of  complex  exponentials  can  be  written  in  the  case  of  the 
gamma  density  for  arbitrary  n. 


B.  Solution  by  a  Matrix  Method.  Any  of  (26a),  (26b)  or  (29)  can  be  solved  by 


o(v. )  =  7-  y  K(y  y.)  o(y. )  ;  i  =  1,  .  . .  ,  N 
i  N  1  j  j 


where  y  =  -1  +  (2k-l)/N.  This  is  a  homogeneous  set  of  N  simultaneous 


(35) 


linear  equations  in  the  p(yj,  i  =  1 . N  and  can  be  solved  by  solving  the 


associated  matrix  eigenvalue  problem 


(K  -  \I)F  =  0 


(36 


Here  P  is  an  N  xl  column  vector  with  components  p(y.),  I  is  the  identity 
matrix,  K  is  the  NXN  matrix  with  element  2K(y.,yJ/N  and  X  denotes  an 
eigenvalue.  The  desired  solution  is  the  one  for  which  X  is  nearest  to  one 
under  the  normalization 


1 


(37 


The  eigenvalue  problem  can  be  solved  by  standard  techniques  from  the  theory 
of  linear  algebra  [18], 


C.  Solution  by  Iteration.  The  integral  equations  lend  themselves  to  solution 
by  iteration  in  several  ways  depending  upon  which  of  the  equations  one  chooses 
to  work  with.  Such  a  choice  will  be  influenced,  in  part,  by  the  kernels  of  the 
equations,  which  are  simplest  in  the  basic  equations  (5a)- (5d)  and  somewhat 
more  complicated  in  (24a),  (24b),  (26a),  (26b)  and  (29).  We  will  describe  in 
detail  only  one  of  these  possible  iteration  schemes  -  that  which  can  be  applied 


to  the  basic  equations. 


With  the  idea  of  iteration  in  mind,  (5a)  and  (5b)  will  be  written  as  1 18] 


Po.n+l(y)  = 


T 

i+y 


d~  pla(">  £01T2"  i*v 


i  - 


(38a) 


Pl.n-l(y> 


l  -  y 


d~  P0nM 


fi  iTZnrr 


(38b) 


r-l 


Here,  p^  (y)  and  p^(y)  represent  initial  starting  estimates  of  p  (y)  and 

p^(y)  respectively.  One  way  of  getting  these  initial  estimates  is  by  using  the 

Fokker-Planck  approximation  p  _(y)  and  (4b)  and  (4c)  in  the  integrands  of 

r  r* 

(26a)  and  (26b);  viz.  , 


poo(y) 


(uQ  +  a, )  T 


alji  f  I  1+- 

<i —  J  ho(T;'’— 


(39a) 


P10<Y)  = 


S^liT 

2  (1  -  y ) 


d"(l-»pFp(- 


1  -  - 
1  -  v 


'-1 


(39b) 


These  then  lead  to  p^(y)  and  p^(y)  by  (38a)  and  (38b),  and  the  results  can 
again  be  used  back  in  the  integrands  of  (38a)  and  (38b)  to  give  p^fy)  and 
p^7(y),  etc.  The  general  procedure  is  illustrated  in  Fig.  1,  which  also  shows 
the  final  step  of  getting  an  estimate  p(y>  of  p(y)  at  the  culmination  of  the 
iteration.  It  is  reasonable  to  suppose  that  (see  [l8]) 


pn(y)  =  Hm  PQn<y> 


n 


(40a) 


P  1  (y )  =  Hm  PlQ(y) 


(40b) 
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and  that  the  convergence  will  be  rapid  if  the  initial  guess  p^^ly)  is  good. 
The  results  of  some  calculations  using  the  iteration  procedure  will  be  given 
in  Sec.  VI.  As  will  later  be  seen,  a  single  iteration  is  sufficient  in  some 
cases. 


Fig.  1.  Illustration  of  iteration  procedure.  The  numbers  above 
the  arrows  are  the  equations  to  be  used  in  going  from 
one  step  to  the  next,  and  p(y)  is  the  final  estimate  of  p(y). 


It  is  also  possible  to  iterate  one  of  (24a),  (24b)  or  (29).  An  example  -will 
later  be  given  in  which  p<y)  will  be  estimated  by  the  single  iteration  of  (2  9);  i.e 


(v>  = 


p '  y 


j  K!y,-lpFp 


r-)d- 


(41a) 


which  is  the  same  as 


p(y)  = 


(U0  +  U1)(1  ^  )J-1 


’  ,  d-  W 


7 ,  (  T  in  — 

ll  1  -  y 


(41b) 


and  it  will  be  seen  that  the  iteration  is  good  to  three  significant  figures.  Such 
accuracy  is  possible  because  the  Fokker  -  Planck  approximation  is  excellent  in 


many  cases. 


VI.  DIFFERENTIAL  EQUATION  APPROACH 

In  certain  cases,  the  integral  equations  (5a)- (5d)  can  be  turned  into 
differential  equations,  and  this  approach  is  the  main  concern  of  this  section. 
First,  however,  we  consider  some  general  relations,  which  can  be  deduced  by 
differentiation,  for  arbitrary  interval  densities. 

A.  Some  Differential  Formulas.  Multiplying  through  by  (1  +  y),  and  differ¬ 
entiating  with  respect  to  y,  (5c)  leads  to 

3uo  17 1(1  pjy}]  =  PqW  "  Pl(^  (42a) 

In  a  Similar  way,  (5 d )  yields 

2  “  1  c^T  t(1  “  y)  P+(y)l  =  P0(y)  '  Pi(y)  (42b' 

Adding  and  making  use  of  (4b)  and  (4c)  then  gives 

3  (u  +  u  ) 

- 2 - -  cTf  ^  '  y2)P(y)]  =  Pofy)  '  /43) 

Equation  (43),  when  combined  with  (12)  leads  directly  to  (17)  and  (22).  This 
derivation  of  (43)  illustrates  the  general  method  of  turning  the  integral  formulas 
into  differential  formulas,  and  will  be  used  in  essentially  the  same  way  in  the 
remainder  of  this  section. 

3.  Intervals  with  Gamma  Densities.  Aside  from  (43),  it  does  not  aopear  pos¬ 
sible  to  use  the  differential  equation  approach  without  specifying  the  interval 
statistics  of  the  input  process,  so  we  here  take  f^(t)  and  f  (t)  to  be  Gamma 
distributed  with  densities  given  by  (32).  Then  it  can  be  shown  by  repeated 
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differentiation  that  (5a)  and  (5b)  lead  to  the  differential  equations 


£+P0(y)  =  (-  a0)n  1  ?1(y) 
p^’y)  =  pQ(y) 


where  Z  and  Z,  are  linear  differential  operators  defined  by 


a  1  -  a 

£  (.  )  =  (1  -  y)  1  (1-y)  (1-y)  (1-v)  1 

dy  dy  dy 


(•  1 


n  times 


3.  1  -  rr 

z  j  •  )  =  (i  +  y)  ~  (i+y)  •••  y) 


(•) 


n  times 


Operating  with  Z  on  (44a)  and  £  x  on  (44b),  it  follows  that 


P0^}  '  ('  aQal)n  1  P0(y)  =  0 
£ pL(y)  -  (-  1  Pj(y>  =  o 


These  are  linear  differential  equations  of  order  2n+2.  When  n  =  0,  we  have 
the  case  of  exponentially  distributed  intervals,  and  it  is  straightforward  to 
show  that  solving  (46a)  and  (46b)  along  with  (43)  leads  to  (9a)  and  (9b).  The 
next  simplest  case  is  that  of  n=  1,  to  which  we  now  turn. 


(44a) 

(44b) 

(45a) 

(45b ) 

(46  a) 

(4bb ) 
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_  _  2  -at  „ 

C.  The  Case  t'(t)  =  a  te  .  For  the  sake  of  simplicity,  we  consider  only 
the  symmetric  case.  Equation  (46a)  for  p^(y)  then  becomes 

1  -  v"  )“  o  4'  ~  2  (2  -  5v~  2a  y)  (1  -  v~  )  o  -  [  (6a“  -  24a  *  25  )v"  -  8  (a- 2  )v  -  (2-“  -  8a  -  5)]p  ’ 

■  0  0  '  ‘  o 


(3-2a)  [(2a“-  6a-5)y-  3  +  2a.]p'  -  [(1  -  a)4  -  a4]pQ  =  0 


(4  7) 


This  equation  is  somewhat  formidable,  and  we  have  not  been  successful  in  find¬ 
ing  general  solutions.  However,  the  form  of  the  a  ■  .-mptotic  behavior  of  p(y) 
around  the  endpoints  v  =  s  1  can  be  determined  from  it.  Setting  x  -  1  -  y, 

1  *  v  =  2  -  x,  d/dy  =  -d/dx  in  (47),  and  assuming  a  solution  of  the  form  x 
leads  to  the  indicial  equation 


v  (v-  l)[v-  (a  +  D]  =  0 


(48) 


Consequently,  a  power  series  solution  to  (47)  will  have  the  leading  terms 

Pq(y)  —  ci  ""  -y)  *  c  (1 -y)  £/i(l-y)  around  y-l  (4Qa) 

Similarly 

Pg(y'~  c_1-c_2(l*y)ic  ^(l+y)a  1  +  c  (1  -y) ~  *  (1  -y)  around  y  =  -  1  (49b) 

The  c's  are  constants  (which  could  depend  upon  a'. 

For  a  =  1,  (47)  can  be  once  integrated  and  leads  to  the  somewhat 
simpler  third  order  equation 


2  2  2  2. 

(l-y  )  Pq  x  2  (2  -  y )  ( 1  -  y  )  p^  x  ( 3  J_y  ^ )  p^  -  (l-v)pn  =  const.,  a  =  1  (50) 


0 


It  is  possible  to  get  a  differential  equation  for  the  unencumbered  density 
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p(y)  in  the  following  way.  Using  (34e)  along  with  (4b)  in  (26a)  leads  to  the 
integral  formula 


p  (y)  -  2a  (l*y)' 


(l--)d- 


-2- 


<5  1 


When  this  is  substituted  back  into  the  right-hand  side  of  (43)  for  p„(y)  and 

J 

p  1  (y )  =  p  (-y),  the  resulting  expression  contains  only  the  unknown  p(y).  The 
integrals  can  be  eliminated  by  differentiation,  which  then  gives  the  third  order 
differential  equation  for  p(y) 

0  =  (l-y2)2p  "-  2(5-2a)y(l-y2)p7  *  2  [(3a2  -  1  1  -  -  1 2 )  y2  -  a2  -  3a-  4  ]  p  ' 

-  2 ( 3 - 2a ) (a 2  -  2a*2)  yp  (52 

As  was  done  above  for  p^(y),  the  form  of  the  asymptotic  behavior  of  p(y)  in 
the  vicinity  of  y  =  =  1  can  be  determined  from  this  differential  equation,  and 
is 

2  a  -  1  2  a  -  1  «  2 

p(y )  ~  k  ^  ~  k^  (1  *  y  )  *  k  ^  (1  ~Y~  )  v/l  d  -  y  )  around  y  =  -  1  1 5  3 

From  this,  it  is  apparent  that  p(y)  -»  ”  as  1  y  !  -*  i  1  if  a  <  1,  and  that  p(y) 
remains  finite  for  all  y  if  a  >  1. 

In  the  special  case  a  =  I  '2,  (45)  can  be  once  integrated,  and  doing  so 
results  is 

[  d  -  )  n  7  ]  -  —  (5v2-3)d  =  constant  (a  r  1  /2 )  (54 

dy  2 

5 

The  constant  can  be  evaluated  in  terms  of  y“  by  integrating  the  differential 


equation  over  [-1,  1],  The  result  of  this  is  the  differential  equation 


—  [  (1  -  v2  r  o  '  ]  -  -7  (5  v"  -  3)°  =  -7  (5y  -  3)  <■ 
dv  ‘  2  ■* 


1  ;->  1 


1  34b 


and  appropriate  boundary  conditions  for  the  solution  are  that  piy>  be  sym¬ 
metric  and  integrate  to  one. 

Just  how  appropriate  boundary  conditions  can  be  found  from  the  integral 
equations,  in  general,  is  not  obvious.  For  the  gamma  dens-tv  with  n  -  - ,  :t  is 
straightforward  to  show  from  the  basic  integral  equations  that  ?  <  1  >  =  phi'll  '  r>< 
But,  the  differential  equation  for  p  (y)  is  of  fourth  order.  For  the  gamma 
density  with  arbitrary  n,  it  can  be  shown  that  p  (y)  and  its  first  n  deriva¬ 
tives  must  vanish  at  n=  1;  however,  since  the  corresponding  differential 
equation  is  of  order  2n~2,  an  additional  n~l  conditions  must  be  speciried 
lone  of  these  is  that  the  density  integrate  to  one). 


VII.  EXAMPLES 

The  results  of  some  calculations  will  be  given  in  this  section  in  the  two 
cases  of  intervals  with  Mcradden  s  density  and  intervals  with  a  gamma  density 
with  n=  1.  Only  symmetric  cases  will  be  considered.  Since  exact  results  are 
known  for  the  McFadden  interval  density,  it  will  be  possible  to  use  them  to  see 
just  how  accurate  the  approximate  techniques  of  Sec.  IV  are.  First,  though,  we 
consider  cases  in  which  the  Fokker -Planck  ’approximation"  is,  in  fact,  exact. 

A.  Me  F  adden  Interval  PDF,  b  =  a  -  I  and  b  =  a  -  2 .  The  McFadden  interval 
density  (1  Oa  /  degenerates  to  the  exponential  density  when  b  =  a  *  1,  and  it  i s 
known  that  the  Fokker- Planck  density  is  then  exact.  ’<Vhat  is  somewhat 
surorising  is  that  the  Fokker  -  Planck  ’’approximation  ’  p„p(y>  is  also  exact 


when  b  =  a~2.  From  (10b),  when  b  =  a  x  2, 


(5  5) 


The  Fokker  -  Planck  result  comes  from  (15),  and  it  can  be  verified  that  (10b)  for 
p(v)  is  the  same  as  (5o)  by  first  changing  variables  of  integration  in  one  of  the 
incomplete  beta  functions  (by  letting  u  ^  1-u)  in  (10b),  and  then  integrating  by 
parts  twice. 


The  probability  density  p(y)  can  be  expressed  as  a  polynomial  in  v  when¬ 
ever  b  =  a  bm  and  m  is  an  integer  -  then  m  integrations  by  part  ere  require 
to  get  the  polynomial  form  from  (10b). 
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and  the  exact  p(y)  is 


«*(t)  =  3e  -  3e  t-  e 


h(t)s3«/Te  ^  sin{VTt) 


,  -1/19 

,  s  _  (1-y2) 

p  (y1  -  ---  — — - 

2  B  (18/1  9,  18/19) 


p<y>  =  ^  (7~y2) 


Table  I  compares  values  of  p(y),  Ppp(y)  and  p(y),  the  result  of  a  single 
iteration  by  means  of  (41b).  As  can  be  seen  from  the  table,  even  though 
pvp(y)  -»  53  as  y  -»  1,  p(y)  approaches  the  correct  value.  Also,  p(y)  is 
accurate  to  three  significant  figures. 


C.  Gamma  Interval  PDF,  n=  1.  For  the  gamma  interval  density  with  n=  1, 
(32),  (33)  and  (16)  give 

2  al+2a  _  ^ 

u  =  7 :  u  =  z  r+  a  :  =  aT 

The  Fokker  -  Planck  approximation  is  especially  simple  when  u=  1;  i.e.,  wh 
a  =  (1  +  a/TT ) / 4 .  Then  Dpp(y)  =  1/2,  and  (39a)  gives 


In  this  case,  (feOe)  and  (lib)  each  result  in  v  =  8  1  /385. 


Table  I.  McFadden  Interval  PDF,  a  =  1,  b  =  4,  ?  =  1. 
Comparison  of  Approximate  and 
Exact  Densities 


y 

p(y) 

p(y) 

PFP  ^ 

0 

.4773 

.4775 

.  4836 

.  2 

.  4800 

.  4801 

.4847 

.4 

.  4882 

.4883 

.4881 

.6 

.  5018 

.  5018 

.  4951 

.  8 

.  5209 

.  5209 

.5103 

.  9 

.  5325 

.5  324 

.  5278 

.  99 

.  5441 

.  5440 

.  5  943 

.  999 

.  5453 

.  5452 

.  6707 

1 

.  5455 

.  5454 
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Poo(y)  =  a 


{(l+y)/Z]ZCL'1  _  ( 1  -ty ) / 2 
(2a- 1 )  (a  -  1 )  a  -  1  2a  -  1 


2  I 


a  =  (1  +  v'l  7 )  /  4 


(62) 


This  p^q(v)  was  used  to  do  one  step  of  the  iteration  as  outlined  in  Fig.  1,  and 
the  results  are  plotted  in  Fig.  2.  Also  shown  on  the  figure  are  the  results  of 
an  8x8  matrix  approximation  to  the  integral  equation  for  p^(y).  As  the  curves 
show,  both  pQ^(y)  and  p^j(y)  agree  well  with  the  matrix  approximation  except 
in  the  vicinity  of  y=  -1.  In  this  region,  only  p^^(y)  has  the  correct  asymp¬ 
totic  behavior  as  predicted  by  (48a).  Also,  p(y)  has  the  correct  asymptotic 
form  as  given  by  (52). 

The  conditional  moment  Afy)  was  computed  from  (17)  using  p^(y)  and 
p(y)  for  pQ  (y)  and  p(y)  respectively,  and  is  shown  in  Fig.  3.  As  the  figure 


shows,  A(y>  is  nearly  linear  except  in  the  vicinity  of  the  endpoints.  By 
(43)  and  the  asymptotic  form  (4aa),  it  can  be  shown  that  p  ^y)/p(y)  -*  4  as 
y  -*  -1,  and  consequently  from  (17)  that 


as  y  -*  1 


(6  3) 


which  is  consistent  with  the  behavior  depicted  in  Fig.  3.  Similar  computations 
were  repeated  for  the  case  a  =  1/2,  and  are  also  shown  in  Fig.  3. 

Fig.  4  contains  curves  of  p  ^  (y)  for  z  =  1/2,  1,  3  and  7  which 

were  obtained  by  the  same  iteration  procedure  as  just  used  for  z  -  (1  ~y/T7)/4. 

These  curves  agreed  to  within  a  few  percent  with  the  results  of  an  8x8  matrix 

approximation  except  in  the  case  a  =  1/2.  Because  p^(y)  tends  to  bunch  up 
around  y  =  -l  for  small  z,  a  matrix  approximation  with  more  points,  or 
with  nonuniform  spacing,  becomes  necessary.  In  contrast,  iteration  was 
found  to  work  well  even  for  the  smaller  values  of  z. 


D.  Me Fadden  Interval  PDF,  a  =  6,  b=4l.  In  all  of  the  previous  cases,  the 
Fokker-Planck  approximations  were  themselves  good  approximations  to  p(y). 

A  case  in  which  this  is  not  so  is  that  of  the  example  used  by  McFadden  in  which 
a  -  d,  b  =  41  and  5=1.  In  this  case, 
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u  =  A  “77  =  1.  99520.  .  .  (64) 

-  n  n  +  6 
n  -  0 
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Fig.  5  gives  a  comparison  of  p(y)  and  p  (y),  and  Fig.  6  shows  the  actual 
and  Fokker-Planck  approximation  to  -A(y)/3^.  The  conditional  moment  is 
far  from  being  approximately  linear,  and  p  (y)  is  not  a  good  approxima¬ 
tion  to  p(y).  Fig.  5  also  shows  some  values  of  p(y)  which  were  computed 
from  (41b)  and  are  listed  in  Table  II.  The  agreement  of  this  single  iteration 
with  the  true  p(y)  is  excellent  except  in  the  region  where  p(y)  is  rapidly 
changing. 

In  doing  the  calculations  to  get  p(y),  it  was  found  easiest  to  evaluate 
the  terms  in  the  integrand  of  (41b)  from 

a  -  L 

3(^  2)  =  1  -  X  (Y)  zk(l-z)b~1‘k  (66) 

k  =  0 

in  which  a  and  b  are  integers,  and 


1  -  3t  „-l  )  B(ST,  2u+l )  sF(s )  / 

2  j  B(sT,  u)  1  -  F(s)  j 


(67) 


where  F(s)  is  given  by  (31).  Eq.  (66)  is  the  result  of  repeated  integration 

by  parts  when  (10a)  is  used  in  the  definition  of  3  (t),  and  (67)  follows  by 

recognizing  that  (39a)  is  a  convolution  when  the  change  of  variable  (7)  is 

3  ^ 

employed  for  y  and  the  change  "  =  -1  +2e  in  the  variable  of  integration. 
The  inversion  in  (67)  was  done  along  the  imaginary  axis  by  means  of  the 
FFT. 
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Table  II.  McFadden  Interval  PDF,  a  =  6,  b  =  41,  3  =  1. 
Comparison  of  Approximate  and  Exact 
Densitie  s 


y 

p(y) 

p(y) 

Ppp^ 

0 

.  5012 

.  5012 

.  6297 

2 

.  5221 

.  5221 

.6177 

.  4 

.  5955 

.  5956 

.  5800 

.  6 

.7206 

.7  297 

.  5101 

.8 

.-*321 

.  3933 

.  3888 

VII.  SUMMARY  AND  CONCLUSIONS 

The  problem  of  calculating  the  probability  density  function  of  the  output 
of  an  RC  filter  driven  by  a  class  of  binary  random  inputs  has  been  studied  in 
detail.  Some  new  integral  equations  were  derived,  and  methods  for  their 
solution  were  developed.  Also,  transformations  of  the  integral  equations  into 
differential  equations  were  investigated.  Exact  and  approximate  results  were 
compared  in  several  examples. 

The  matrix  solution  technique  was  seen  to  work  well  for  the  cases  in 
which  the  densities  being  computed  were  smoothly  varying,  but  requires  a 
finer  grid  with  more  points  to  yield  any  erratic  behavior.  Iteration  was  seen 
to  work  well,  with  only  a  single  iteration  required  when  the  initial  estimate 
was  good.  In  the  last  example  of  the  McFadden  interval  PDF  with  a  =  6, 
b  =  41,  even  though  the  initial  estimate  was  not  particularly  good,  the  single 
iteration  was  within  10%  of  the  correct  value.  This  example  and  the  first 
example  are  felt  to  establish  the  veracity  of  the  iteration  process.  The 
Fokker  -  Planck  approximation  was  seen  to  not  only  be  a  good  estimate  for 


starting  the  iteration,  but,  in  many  cases,  to  be  in  itself  an  excellent 


approximation  to  the  desired  p(y). 

All  in  all,  the  differential  equation  approach  was  disappointing  since 
it  had  been  hoped  that  it  would  lead  to  new  closed  form  solutions.  The  differ¬ 
ential  equations  are  not  only  complicated  in  themselves,  but  correct  boundary 
conditions  for  their  solutions  are  not  apparent  in  the  general  case.  Never¬ 
theless  some  useful  asymptotic  information  was  obtained. 


-0.5 


0 
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The  density  pQ^(y)  for  a  gamma  interval 
PDF  with  n  =  1  and  1  as  a  parameter. 
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Output  probability  density  functions  for  McFadden  interval 
PDF  with  a  =  6,  b  =  41  and  3=1. 
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McFoaoen  nfervoi  =CF 

C7  =6,  0  •  41 


Exact  and  approximate  conditional  moments  for  a 
McFadden  interval  PDF  with  a  =  6 ,  b  =  41  and  3  =  1. 


APPENDIX.  CALCULATION  OF  yc 


This  appendix  shows  how  y“  is  determined  from  the  system  equation 
for  arbitrary  interval  statistics  of  the  binary  input  process.  We  will  give  the 
derivation  only  for  the  symmetric  case,  and  rely  heavily  upon  some  results  in 
Cox's  book  [3]  on  renewal  theory. 

The  mean  square  of  y(t)  can  be  obtained  by  integrating  the  power 
spectral  density  S^(f).  Relating  S^_(f)  to  the  spectral  density  S  (f)  of  the 
input  binary  process  then  results  in 


2 

y 


d  r 

2  - 


(A 


S  (f)  can  be  further  expressed  in  terms  of  the  generating  function  of  the 
number  of  renewals  of  an  equilibrium  renewal  process  as 


S  (t) 
x 


R  (rHe1'"*  *  e  ” 1 ''  ~  >  d  ’ 
x 


G  (-iai  -1)  +  G  (i  r  -  1 ) 

e  ’  e  ’ 


(a; 


where  R^  (")  is  the  autocorrelation  function  corresponding  to  (f\  and 
G  (s,  \)  is  the  generating  function  defined  by  Cox  and  evaluated  by  him  as 
[3,  eq.  (3.2.6)] 


G  (s, 
e 


s 


2  1  -  Ffs) 

2  1  -  F(s  1 


( A3 ) 
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